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ABSTRACT 

Microlensing light curves are typically computed either by ray-shooting maps or 
by contour integration via Green's theorem. We present an improved version of the 
second method that includes a parabolic correction in Green's line integral. In addition, 
we present an accurate analytical estimate of the residual errors, which allows the 
implementation of an optimal strategy for the contour sampling. Finally, we give a 
prescription for dealing with limb-darkened sources reaching arbitrary accuracy. These 
optimizations lead to a substantial speed-up of contour integration codes along with a 
full mastery of the errors. 



Subject headings: gravitational lensing - methods: numerical - binaries: general 
etary systems 
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1. Introduction 

Mi crolensing is one of the mo st promising methods for finding the first Earth-like extrasolar 
planet (lGaudill201Cll : lDominikll2010l ). When a compact object transits very close to the line of sight 
of a background source star, the flux coming from the source is amp lified by gravitat ional lensing 
and follows a typical bell-shape light curve, analytically described by iPaczyiiskil (Il986l ). If the lens 
is a star accom panied by a secondary b ody like a planet, an additional bump or dip appears on 



the light curve (jMao &: Paczviiskilll99ll ). The timescale of such featu res ranges from a few hours 
to a few days, depending on the square root of the mass of the planet (jGould &: Loeblll992l ). Such 
short timescales require intensive monitoring by telescopes situated all over the Earth. 

Nowadays, more than 30 telescopes are involved in microlensing searches towards the Galactic 
bulge. More than 600 events are discovered every year. Out of these, roughly from 10 to 20 events 
show anomalies that can be interpreted as due to binary lenses. Some of these are finally accepted 
as showing evidence of an extrasolar planetary system. Since the start of microlensing searches. 
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26 events have been reported as containing planetary candidates with stronger or weaker evidence 
DominikI (j2010l ) . Nine of these events have a lso been included in the most upda.ted exoplanet list 



available on the web (http://exopla net.eu) (I Bealieu et al 



2004 



Udalski et al 



Dong et al.ll2009l : iGaudi et al.ll2008. : Gould et alJ 12009 : Ijanczak et al.ll2009l : ISumi et alJl20ld : 



2006 



Bennett et al 



2004 



Bond et al 



20051). Final ly, one event brings the spectacular signature of two planets in the same 



system (jGaudi et al.l 120081 ) . 



Unfortunately, interpreting binary microlensing events is a very long process that may even 
take from one to four years for a single event. This is due to several reasons: one is related to the 
difficulty of getting rid of all systematic errors in the photometry of each dataset. Some sets of 
images need to be reduced with different methods in order to compare the effects of systematics, 
and subtract them from the final result. After the reduction process is complete, the final datasets 
may contain hundreds or thousands of data points, which are then ready for the modelling process. 

The modelling process is typically driven by the criterium of minimization, which can be 
achieved either by downhill algorithms or by Markov Chain Monte Carlo (MCMC) methods. In 
order to evaluate the \^ for a single tentative model, it is necessary to compute the microlensing 
magnification at each data point. However, as is well known, microlensing magnification cannot 
be calculated analytically if the lens is a binary system. As the angular extension of the source 
plays a major role in the observed magnification, the computation of even one single model point 
is relatively time-consuming. Multiplying this basic time unit for the number of data points, the 
number of models within a MCMC simulation, the number of different hypotheses to be checked 
(parallax, xallarap, orbital motion, limb darkening, binary source, . . .), we can easily imagine why 
the study of a single event takes so long. It is then critical to reduce the computational time of a 
single model point as much as possible, so that the whole modelling process is cut down to a more 
reasonable duration. 

There are basically two classes of methods for the computat ion of the microlen sed flux of a 
source. The first is based on the construction of ray-shooting maps (iKayser et al.lll986l ). In practice, 
rays are shot back from the observer to the lens plane and then deflected to the source plane. If 
they intercept the source disk they are counted as contributing to the total magnification. This 
method has three main advantages: it is conceptually simple, it can naturally take into account the 
limb darkening profile of the source, maps at fixed lens configurations can be re-used for different 
source p ositions. Numerous optimized versions have appeared in the li terature, designed for re-use 



of maps (jWambsgansd 119971 : iRattenburv et al 



2002 



Dong et al.ll2006l ). An alter native strategy is 



to sp eed-up the computation for a single source position giving up the map re-use (jBennett &: Rhie 
19961 ). In this case, starting from the positions of the centers of the images, one shoots rays only 
where really needed. A further improvement of this method employing a polar coor dinate-gr i d wit h 
an optimized prescription to handle limb darkening has been recently presented by lBennettl (|2009l ). 



The second method is based on an application of Green's theorem (which can be viewed as 
a two-dimensional version of Stokes' theorem). In practice, one can find the area of an image 
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by calculating a Rie mann integral along the ima ge contour. The use of contours in gravitational 
lensing d ates back to ISchramm &: Kavseii (|l987l ). Green's t heorem w a s the n used for calculating 
areas bv lDominikI (1199311. The method h as be en refined by iDominikl (|l995l ) and then applied to 
microlensing bv lGould &: Gaucherell (|l997l ) and iDominikI (|l998l ). The appeal of this method is that 
a two-dimensional calculation is turned into a one-dimensional calcul ation, which i s mu ch faster, 
in principle. Related to this approach ar e the algorithms presented by iDong et al.l (j2006l ) and the 
adaptive grid search by iDominikI ()2007l ). The main advantages come from the potentially high 
computational velocity and the high flexibility for models in which the lens configuration changes 
(e.g. in the treatment of planet orbital motion). However, the contour integration approach re- 
quires an images reconstruction procedure (which can be sometimes complicated); in addition, limb 
darkening cannot be naturally incorporated in the algorithm. In particular, the latter limitation 
has oriented the community to give a general preference to ray-shooting methods. 

Nevertheless, apart from its undisputed elegance, the contour integration approach is still 
competitive for obtaining preliminary microlensing models very quickly, which is particularly inter- 
esting in view of the realization of real-time modelling of binary microlensing events. Furthermore, 
when orbital motion is relevant, traditiona l ray-shooting rnethods typically bec ome definitely too 
heavy. In this case, only adaptive methods (jBennett &: Rhidll996l : lBennettll2009l ) can compete with 
Green's theorem algorithms. 

Finally, besides ray-shooting and contour integration methods, it is worth mentioning that 

when the source size is only marginally relevant (for sources not to o close to caustics), one can ap- 

proxi mate its effects by quadrupole or hexadecapole approximations (|Gouldll2nn8l : IPeicha O. fc Hevrovskvl 
20091 ). These methods allow to obtain a substantial speed-up of the code avoiding useless heavy 
computations when the source size correction is small. They can be used in combination with other 
methods that may intervene when the source gets closer to a caustic. 

In this paper we present four new ideas for boosting codes based on contour integration ap- 
proach. In Section 2 we show how Green's line integral can be approximated to third order in- 
troducing a parabolic correction, with a substantial improvement in accuracy. In Section 3 we 
present accurate estimates of the residual errors in Green's integral. These are used to implement 
an optimal sampling strategy that allows to minimize the calculations for a given required accuracy, 
as explained in Section 4. In Section 5 we suggest an easy prescription for the treatment of limb 
darkening that achieves a fixed accuracy avoiding lengthy calculations. The benefits achieved by 
all these innovations are documented by several numerical examples in Section 6. 
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2. Green's line integral to third order 

2.1. The concept of Green's line integral 

Consider a generic continuous gravitational lens mapping between the image plane x and the 
source plane y 

y = m. (1) 

Consider a circular source As with radius centered in the position ys- The boundary of the 
source is a circle of radius p* that we shall indicate by 75. A trivial parametrization of this curve 
is 

^//iN ^ I cos 6 \ 

\ smb* / 

For each 9, we can solve the lens equation ([1]) with y = y{6). As 9 runs from to 27r, the 
solutions of this equation describe several curves 7/ in the image plane, parameterized as xj{9). 
The subscript / runs from 1 to the number of images A^. In the case of a binary lens, A = 3 if 
the source is outside all caustics, A^ = 5 if the source is completely inside a caustic. If part of the 
source is inside a caustic, then two images are created at some 9c and disappear at some 9d, so 
that N = 5 with two images xi{9) defined only in the subinterval [9c, 9d]- Creation-destruction 
of images may also occur in several disjoint subintervals of [0,27r], if the source touches two or 
more caustics. All curves 7/ have definite parity = ±1 and represent the boundaries of the 
regions in the image plane that are mapped to the source As through the lens map /. Such regions 
represent the physical images of our source. The ratio between the total area A of all images and 
^4^ represents the sought magnification factor. 

By Green's theorem, the area enclosed by a closed curve 7 is 

A = ±^ j xAdx, (3) 
7 

where the positive sign is taken for counterclockwise curves and the negative sign is taken for 
clockwise curves. We remind that the wedge product between two vectors is a pseudoscalar in two 
dimensions: x Ay = xiy2 — X2yi- 

As our parametrization of the source boundary 75 is counterclockwise, positive parity 7/'s are 
still counterclockwise, whereas negative parity 7/'s are clockwise. Therefore, the total area of all 
images can be found as 

A = ^ -p/ / XI A dxi. (4) 
/ 

7/ 



Such expression still holds also when part of the source is inside a caustic ()Dominikl Il995l ; 
Gould Gaucherellll997l l. 



- 5 - 



2.2. Trapezium approximation of Green's integral 

In order to find a numerical approximation to Eq. we must introduce a sampling of the 
source boundary in the following way 



ys + P*\ . ^ I , (^) 




where {9i} is an arbitrary ordered sequence of n numbers with = < < ■ ■ ■ < < • • • < = 
27r. One simple possibility is to take a uniform sampling — 9i = const. However, this is not 
necessary and more optimal choices are possible, as will be explained in Section SI 

For each 9i, we solve the lens equation ([T]) and find the corresponding points xj^i on the image 
boundaries 7/. If is close enough to xj^i, it makes sense to approximate Eq. (jU as 

^ n— 1 

I i=0 



I i=0 
^ n—1 

- ^ {Xl,i+1,2 + Xl^i^2) {Xl,i,l - , (6) 

where the last version is simply the trapezium approximation of the Riemann integral of the function 
xifi{xi^i). It is more advantageous numerically in that it has one multiplication instead of two. 

Eq. dll) is written in the case of no caustic crossing. It can be easily extended to the general 
case by letting i run only on the values for which the image 7/ exists and adding up connection 
terms between each pair of created images and each pair of destroyed images (see also Section 12. 5p . 

Summing up, in the implementation of the trapezium approximation of Green's line integral 
([6]), we need the following routines: 

• A routine solving the lens equation f or each source pos ition yi. For example, one can use the 



zroots routine of Numerical Recipes (jPress et al.l 120071 ) 



A routine associating the solutions xj^i+i found at each {i + 1)— th step with the correct 
image 7/. This can be done by re-ordering the solutions in such a way that \xi^i — < 
\xi,i — xj^ij^i\ for each J ^ I. If new images are created or destroyed, they will be recognized as 
the last two unmatched images. Of course, this association routine has some failure probability 
when the new solutions are too far from the old ones xi^i. However, we will see in Section 

[3] that a careful estimation of the errors will easily recognize such situations. 



We find that roughly 80% of the machine time is spent in the root finding routine, for which 
there is basically no hope of further optimization (we already re-use old roots as starting values 
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for the next calculation). So, the only way to speed up a contour integration code is to reduce 
the number of points in the sampling while keeping the same accuracy. This can be achieved by 
pushing the numerical approximation of Green's integral to higher orders. 

Another poss ibility to get a round the problem of root finding is the use of adaptive grids 
on the lens plane iDominild (120071 ). In these algorithms, the sampling of the image boundaries is 
obtained by a grid construction directly on the lens plane. Although we will mostly refer to the 
scheme described in this subsection (source sampling and lens equation solving to obtain an image 
sampling), most of the concepts introduced in this paper can also be applied to algorithms based 
on direct sampling on the lens plane. 



2.3. Parabolic correction of Green's integral 

Going back to Eq. Q, we can write it as 



27r 



A = Y,Ipi [ ^ x'ldO, (7) 
^ 

where the prime denotes derivation with respect to the parameter 9. 

Let us consider the generic image 7/ and the generic interval with size A9. The 

contribution of this interval to the whole integral is 

ei+A0 

dAi = ^ j xiAx'jde. (8) 

e. 

The first order trapezium approximation used up to now just reads (see Eq. ([6])) 

dAf'> = hii9i)Axi{ei + Ae). (9) 

Comparing the expansions of dAj and dA^p in powers of A^, we find that they coincide at the first 
and second order, the difference being of order A^'^. 

Now, let us introduce the following correction term 



dA 



ip) 



1 



^/ - ^ L '^'^ ^ 1^. + ^ I^»+aJ (10) 

Adding this correction to the trapezium approximation and comparing the power expansion to that 
of the exact integral ([8]) , we have 

dAi = dAf + dA^ + 0{Ae^). (11) 



The residual error is now of order A9^, which is much smaller than what can be achieved by 

(p) 

the trapezium approximation. dAj can be viewed as a parabolic correction as it takes into account 
the local curvature of 7/ stored in the second derivative x". 
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2.4. Implementation of the parabolic correction 



dAY^ is expressed in terms of derivatives witli respect to 6 calculated at 6i and In 
principle, these derivatives can be easily calculated an alytically in terms of local quantities. The 
explanation is easier if we switch to complex notations (|Wittlll99Cll ). The coordinates in the source 
and lens planes respectively become 



C = yi + iy2 

Z = Xl + iX2- 



(12) 
(13) 



The lens equation for a binary lens with mass ratio q and separation a assumes the form 

1 



from which we obtain 



l + q\z + a/2 z-a/2j ' 



(14) 



dC/dz = 1 
dC 1 



1 



+ 



dz l + q\{z + a/2f (z - a/2)2 
2/1 



+ 



l + q\{z + a/2f {z-a/2) 



(15) 
(16) 

(17) 



The Jacobian determinant is 



J = 1 



dz 



(18) 



Note that the Jacobian determinant must be calculated in the linear approximation too, in order 
to assess the parity of the image. 



Deriving Eq. (|14|) with respect to 0, we have 



OZ 



(19) 



Inverting this equation with its complex conjugate, we get 



oz 



J- 



(20) 



Deriving again with respect to 0, we get the expression for z" 

32 



z 



^ dz^ ^ ' dz 



J 



-1 



(21) 
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The key fact is that ah these quantities can be calculated exactly starting from the source 
parametrization, which in complex notations reads 



C = C5 + P*e^^ (22) 

where Cs = ys,i + Ws,2 is the center of the source disc. From this expression, we get = ip*e.^^ 
and C" = -P*e'^. 



The parabolic correction (jlOp contains terms of the type 

x' Ax" = - iz"z' - z'z") . (23) 

Plugging the former expressions for z' and z" and using the source parametrization, we finally 
obtain the compact expression 



x' A x" = s + Im 



^ ^ dz^ 



J-\ (24) 



The implementation of a parabolic correction is therefore relatively simple. For each ^, after 
the extraction of the roots of the lens equation, we just have to calculate z' , J and d^jdz^ for 
each root, taking = ip^^e^^ , and then store the value of the wedge product (f24l) . 



Finally, when we compute Green's line integral, for each arc [x/^j, x/^j+i] we can put together 
all the ingredients to calculate both the trapezium approximation ([9]) and the parabolic correction 

(HOD. 



2.5. Parabolic correction at critical points 



As pointed out before, it might happen that a portion of the source boundary lies inside a 
caustic. In this case, at some 6c a pair of new images is created and at 6d another pair is destroyed. 
Green's theorem can still be applied, but since 6c and 9^ do not generally belong to our sampling 
{6i}, we need to introduce connection terms between the starting points of the created images 
(the same happens for the pair of destroyed images). Let us discuss the case for pair creation, the 
destruction being analogue. 

For a given sampling {^i}, the new pair of images appears at some 9i, with 6i-\ < 6c < 6i. Let 
us call the starting points of the new images x+^j and x_,i. The problem is that the parametric 
distance between the two images is not available, since the precise value of 6c is unknown. However, 
the standard expansion of the lens equation in a neighborhood of a fo ld tells us that the two created 
images move away from the creation point as ([Schneider et al.lll992l ) 



x± 



X 



-byi±^2ay2X2+{b^-ac)yf 



(25) 
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where A, a, b and c are coefficients of the expansion of the lens equation near a fold, (^1,^2) is the 
position of the source relative to the fold caustic and x± is the position of each of the two images 
relative to the critical curve. 

Expanding our parametrization in the neighborhood of the crossing point, we have y = {6 — 
^c)(ci, C2), with ci and C2 being two constants depending on and 9c- The sought connection term 
reads 



dAc 



J x+ A x'j^dt - J X- A x'_dt, (26) 



where we have assumed that x+ is the positive parity solution while x_ is the negative parity one. 
The trapezium approximation with the correct signs is simply 

= ^ {X-A,2 + X+,i,2) {X-,L1 - . (27) 

Now, we propose the parabolic correction 

^4^^ = ^ [(^V. A f'l,) - (x'_, A x'L,)] m\ (28) 

where 

Ke = ' ^+'^ ~ (29) 



replaces the parametric distance used in the ordinary parabolic correction. 

Using the approximate general expressions for the images (j25p . and expanding dAc, dAc^ and 



dA^^ in powers of {9i — 6c), we realize that the trapezium approximation is accurate only to first 
order in {9i — 9c), the residual error being of order {9i — 9c)^^'^- The parabolic correction accounts for 
the term of order {6i — 9c)^^^ and leaves a residual error of order {9i — 9c)^l'^ . Note that the orders 
of the errors in arcs containing a pair creation or destruction are halved with respect to ordinary 
arcs. This is a major reason for increasing the sampling of the source boundary near caustic points. 
With a uniform sampling, instead, the error would be largely dominated by intervals containing 
caustic crossings. 



2.6. Final remarks on the parabolic correction 

The parabolic correction has some computational cost because it requires some additional 
operations to be performed on each root. However, such cost remains negligible with respect to the 
time spent in the root inversion routine. Moreover, it helps reducing the residual error dramatically 
for each sampling interval up to the fifth order in A0. We can read this achievement in 

the other way round: with the parabolic correction we can reach the same accuracy as with the 
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trapezium approximation only, but with a much sparse sampling of the boundary curves. Since, 
as said before, most of the computational time is spent in the root inversion routine, which must 
be run once for each 6i, we can aim at a substantial speed up of the code if we are able to reduce 
the sampling as much as we can without losing accuracy thanks to the parabolic correction. A 
fundamental step toward this goal is a careful estimation of the residual errors, which is the subject 
of the next section. 



As shown in the previous section, the residual errors in Green's line integral after the introduc- 
tion of the parabolic correction are of the order A^^. The exact expression for the fifth order term in 
the power expansion of Green's integral ([8]) contains third derivatives of z, as can be easily guessed. 
However, we do not want to make more calculations for the estimate of the fifth order term and 
rather use the quantities already calculated to make a realistic but economic estimate. Secondly, 
higher and higher order derivatives are more and more affected by numerical errors. Finally, we 
must also take into account the possibility of a wrong matching of the images, as anticipated in the 
previous section. 

For all these reasons, we disregard the fifth order term in Green's integral and prefer to 
introduce three new quantities as error estimators. These three quantities are built up from first 
and second derivatives of z, thus requiring a minimum amount of additional calculations. They are 
meant to intervene in different situations with the aim of being complementary to each other and 
cover all possible sources of error. 

The first estimator is 



As the parabolic correction is based on an average of the wedge product x'j A x'j- on the two 
end points of the arc, it is natural to estimate the error using the difference of the two quantities 
that are averaged. A power expansion in AO reveals that Ei is of order rather than A9^. 
The main reason for using -E/ j i is that it performs very well in the identification of wrong images 
matching, since the two wedge products take very different values if the two points do not belong 
to the same image. On the other hand, it does not seem to systematically dominate over fifth order 
error estimators. 

The second estimator is 



3. 



Error Estimators 



3.1. Error estimators for ordinary images 




(30) 




(31) 
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The ratio of the squared distance between the two points and the scalar product of their derivatives 
is an approximation to (see the definition of in equation ()29p ). It can be shown that -£7,2,2 
is of order It also takes large values in case of wrong images matching but works in a 

complementary way to Ei as it is based on different quantities. In particular, this estimator is 
particularly effective for the detection of hidden cusp crossings between yj and y^+i, which might 
otherwise be a very dangerous situation. 

The last estimator is 

A^^ (32) 

which intervenes in undersampled situations when -E7,i,2 is incidentally zero. 

The error estimators just defined can be combined so as to build an error estimate for each 
arc [x7,j, X7,i+i]. We adopt a simple sum 

Ei^i = + -E'/,i,2 + -E'/,i,3- (33) 

We prefer a simple sum to the usual quadrature combination of the errors in order to minimize the 
operations while keeping a more conservative attitude. 



dAf 



3.2. Error estimators at critical points 

For the connection terms between image pairs created at some critical points we need different 
error estimators. 

We define 

^i'^ = 4^ I (^'+. A f;,) + (x'_,, A x'L,) I a/, (34) 

in analogy to This quantity is of order 2 in {6i — 9c), instead of order 5/2. The same 

comments as for Ei apply. 

The second estimator is 



AO, (35) 



which is still of second order in (9i — 9c)- The upper sign applies at creation of two images and the 
lower sign applies for destruction. Indeed, at the creation of two images the two starting points are 
expected to move far apart in opposite directions. Conversely, at destruction of two images, the 
end points converge to the same critical point, e!^-' becomes very large if this is not the case. 



The third estimator is 

10 

which is analogous to -E^/,i,3 and is of order 5/2 in {9i — 9, 



AO^, (36) 
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3.3. Resurrecting buried pairs of images 

Finally, a very dangerous situation may occur when working with Green's theorem approach 
and a very sparse sampling, which might be likely if we take full advantage of the parabolic cor- 
rection. If the source boundary grazes a caustic, it might happen that a very small slice of the 
source is inside it. This slice might entirely fall in the middle between two sampling points, so 
that 9i < 6c < Od < Oi^i. In this situation, we have no idea that the i-th. interval contains an 
additional pair of images. These images could be completely missed with a consequent dangerously 
large error, whereas all remaining images are calculated to the desired precision. 

In order to thwart this menace, we introduce an additional estimator in the following way. 
Suppose that at 9i the source boundary is outside the caustic. Then we have three real images 
satisfying the complex lens equation plus two additional roots that satisfy the fifth order polynomial 
version of the lens equation but not the original one. When the source crosses a caustic, these two 
ghost roots merge into a double root and then become real. 

Therefore, we can estimate how far yj is from a caustic evaluating the distance between these 
two ghost roots 

9i = \Zgl,i - Zg2,i\. (37) 

Our idea is to check whether gi^i — gi> gi- By linear interpolation one would expect that at step 
i + 1 the source is inside the caustic and the two ghost roots have become real. If this is not the 
case, we add an error 

EG,^ = {g^-l-9i?, (38) 

which might be quite large. As we shall sec in the next section, this large error will drive the 
optimal sampling strategy to oversample the i-th interval in search for a possible caustic crossing. 

Of course, for symmetry, we also check that gij^i —gi > gi and add an error Ec^i-i = (ffi+i —gij^ 
if the source is outside the caustic at step i — 1. This check on the ghost images turns out to be 
very effective in discovering buried caustic crossings in undersampled regions. 

4. Optimal sampling 

A uniform sampling of the images amounts to adopting a fixed step size AO, so that 0^+1 — 9i = 
AO for each i. A uniform sampling may be quite inefficient: in fact, it might oversample regions 
with little contribution to the magnification or with small errors whereas caustic crossings, which 
require a denser sampling, would not receive any particular regard. 

For these reasons, we adopt a different sampling strategy. Consider a given n-points sampling 
{^1, . . . , On}, with ^1 = and On = 27r. The question is where to put the next sampling point 6. 



- 13 - 



For each interval ^i+i], an estimate of the errors is given by 

Ei = Y,Ei,u (39) 

where I runs on all images present in the i-th interval and the error estimators Ej^i are defined 
in Section [3l If any images are created or destroyed in this interval, we also add the errors of 
the corresponding connection term. Moreover, we also perform the additional check described in 
Section [3.31 and add the corresponding error term if the source is grazing a caustic. 

Then, we select the interval (labelled by i) with the largest error (£"• > Ei for all i ^ i). 
Finally, we simply set the new sample point as the midpoint of the interval with the largest error: 

The sampling is thus increased only where really needed, starting with the intervals with the 
largest errors. Since the new sample point lies in the middle of the sampling sequence, once we 
calculate the images x{9) of y{6)^ we need to make the correct association both with the images 
x{9{) that precede x{9) and with the images x(^_,_i) following, so as to have the new re-sampled 
boundary curves. After the association is done, we can recalculate the magnification contributions 
and the errors in the new sub-intervals [0i,9] and [6, 0~^j^^. As can be easily guessed, this procedure 
is technically easier to achieve by defining the boundary curves and the sampling sequence as linked 
lists rather than arrays, so that we can easily cut and link them so as to insert new members in 
the middle. 

As a starting minimal sequence, we take {0, vr, 27r} with only two points {9 = 27r being just a 
replica 9 = 0). The third sample point will thus be 7r/2 or 37r/2 depending on the respective 
errors of the two initial arcs, and so on. Iterating our sampling procedure, we will end up with 
a non-uniform sampling, with more points where the errors tend to stay larger (typically close to 
caustic crossing points). If we are far from any caustics, the sampling will tend to be uniform, 
anyway. 

The full control of the errors allows us to establish when to stop the iteration. In fact, the total 
error in the magnification is simply given by the sum of the errors of all intervals in the sampling 

E = Y,Ei. (40) 

i 

Therefore, if our target accuracy is we just have to iterate until 

E/{T,pl) < (41) 

In some situations (at high magnification points), the images are very thin stretched arcs, in 
which the errors of the inner side compensate the errors of the outer side. In these cases, thanks to 
this cancellation, the true error is much smaller than the sum of the absolute values of the errors 
of all intervals. Our algorithm can thus be stopped earlier than prescribed by Eq. (j4ip . More 
precisely, we stop when the magnification has changed by less than 5ijl/2 in the last n/2 sampling 
steps. 
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5. Limb Darkening treatment 



The main drawback of codes based on contour integration approach is that the source is treated 
as a uniform brightness disk. Real stars are actuahy non-unifor m, with a brightness profile that 



can be approximated by a linear limb-darkening law (jMilne 



1921 



I{p) = Ifip/p.) 

m 



1 - a ( 1 - a/T 



1 - a/3 L 



where r = p/p^ and / is the average surface brightness, so that 

P» 1 
j 2ti pl{p)dp = 1 j 2rf{r)dr = I. 
* 

A typical profile is shown by the dotted line in Fig. [TJ 



(42) 



(43) 



Indeed, whenever a caustic crossing is present in microlensing, limb darkening cannot be ne- 
glected in accurate modelling of the event. Therefore, if we want to use our Green's integral 
approach with the parabolic correction in real microlensing events, we must find a way to incor- 
porate limb darkening. One obvious solution is to calculate the magnification of several concentric 
disks at the source position ys with radii given by {pi, • • • ,Pm}) with Pm = P*- Each disk should 
be weighted according to the limb darkening profile in order to build up an approximation to the 
correct result. In this scheme, for a single point in the light curve, we must calculate the magni- 
fication of m disks instead of just one. As a consequence, the computation is slowed down by a 
factor of m. 



An alternative possibility, proposed by iDominikI (jl998l ) is to implement Green's theorem with 
different integrand functions taking limb darkening into account. In any case, additional integra- 
tions are required from the center of the source to the periphery. 

Unfortunately, there is no known way around this problem. Our approach follows the standard 
multi-disk solution improved by a clever choice of the radii pi and by a full control of the errors. 



5.1. Magnification and errors in annuli 

Let us consider an annulus of our source with inner radius pi-i and outer radius pi. The total 
luminosity of the source annulus is given by 

Pi 

lf)= I 27rpI{p)dp = l7rpl[F{n)-F{u.i)], (44) 

Pi-i 
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where ri = pi/ and we have introduced the cumulative function 

r 

; J dr'r'f{r'). (45) 



Fir) = 2 

b 



Gravitational lensing introduces a point-source magnification factor p{r^ 6) which modifies the 
observed luminosity as 



27r 



Ii = Ipl j rf{r)dr j d9p{r,9). (46) 



Ti-l 



The total observed luminosity is given by the sum of the luminosities of all annuli. Dividing 
by the original source luminosity Iirpl, we get the limb-darkened magnification factor 



^ m m 
i=l i=l 

Mi = ^ j rf{r)dr j d0/x(r,0). (48) 
n-i 

Using our contour integration approach, we are able to estimate the magnification factor for a 
uniform disk of radius pi to an arbitrary accuracy Sp,. Of course, such a magnification factor for a 
finite size source is just the average of the point-source magnification p{r, 0) on the source disk 

l^i = :^ j rdr j dep{r,e). (49) 



Eq. (j49p looks very similar to Eq. (|48p . save for the profile function /(r) appearing inside the 
radial integration in Eq. (j48p . Therefore, we can approximate the contribution Mi of each annulus 
to the total magnification by replacing the brightness profile /(r) by a constant average brightness 

F(ri) - F(ri^i) , , 

fi = (50) 

Taking fi out of the integral, we get the following approximation for Mj 

Mi = fi [pirf - p^^lrtl\ (51) 

in which all objects involved can be easily calculated in our code. The approximate expression Mj 
reduces to the exact one Mi for very thin annuli. As we can see in Fig. [H the linear limb darkening 
profile is approximated by a block function, in which each block has a constant brightness given 
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by fi. Increasing the number of bins, the hmb darkening profile is approximated better and better. 
More specifically, if 5r = ri — rj_i, the difference between Mj and Mj is of the third order 



2-K 



5Mi = -^f'{n)5r^ / dedrii{r,e). 



(52) 



This expression for the residual error can be used to construct an efficient and economic error 
estimator without calculating derivatives explicitly, namely 



(1) 



1 



(53) 



Note that 5M^^^ reduces to 5Mi only if we neglect the second and higher derivatives of //(r, 9) 
in a neighborhood of the source. This means that this error estimator could be unreliable in some 
situations in which /i(r, 6) has a high curvature. We will come back to this issue later. 

At caustic crossings, /x(r, 9) diverges and 5Mi loses meaning. In principle, our error estimator 
6M^^^ does not diverge but does not track the error correctly. Therefore we introduce a new 
estimator that should be used whenever the number of image contours at rj_i differs from the 
number of contours at 



6M„ 



1 



[rffJ-i - rf_ifJ.i-i] [/(n) - /(ri_i)] 



(54) 



which is always regular and of order 6r in the limit of thin annuli 



5.2. Sampling the source profile 

We have now an approximate form for the magnification of the annuli and error estimators to 
control the accuracy. We must now give a prescription for the choice of the radii of our annuli in 
order to complete the limb darkening treatment. 

Starting from a sequence ri, • • • , r^, we select the annulus with the largest error, say [r7_i, r^]. 
Then we divide it in two annuli, by inserting the radius f between 'ri_i and rj. The new radius f 
is chosen in such a way that F^t^) — F{f) = F{f) — F(r^„]^). In this way we make an equipartition 
of the cumulative function. As a practical example of this partition criterium of the source, in Fig. 
[T]we show a linear limb darkening profile together with a block approximation with four bins and 
a block approximation with 16 bins. The radii are chosen so as to have F{ri) — F(rj_i) = l/riiyins 
and the constant brightness value in each bin is given by Eq. (j50p . We can see that this block 
approximation rapidly converges to the exact profile when the number of bins is increased. Thanks 
to our error control, however, we do not need to increase the sampling everywhere but only where 
really needed. For example, if only the periphery of the source intercepts a caustic, our error 
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estimators will require more annuli to be created close to p* without calculating useless annuli at 
the center of the source. 

Finally, let us come back to the issue of the second derivative of fi{r,6). It might happen 
that the finite size magnifications of the disks steadily grows from the center to the periphery. 
However, it might also happen that at some radius rj Hi starts to decrease. In such a situation, 
it might happen that — /ij leading to a dangerously small 6Mj^\ In order to overcome this 
problem, when adding a new radius f to the sampling between rj_^ and r^, with its finite size 
magnification p,, we calculate the errors of the annuli [r^_i,f] and [f,rj\ according to Eq. (|53p or 
(j55p and then add to both annuli an error 



which accounts for possible changes of slope in the finite size magnifications /ij. 

Summing up, with this error-driven sampling strategy we continue adding annuli until the 
total estimated error drops below the desired accuracy (5/i. Each finite size magnification fii is also 
calculated at accuracy 6fi. Since each of them is weighted by the average flux in the expression of 
Mi, the total error in M coming from the /ij is kept below 6fi. 



In this section we will consider some explicit examples of magnification computations with the 
aim of illustrating the power of the innovations proposed in the previous sections. 



In order to present a test as exhaustive as possible, we calculate magnification maps with 
different levels of target accuracy (5// and evaluate the relative difference. We take maps calculated 
at 6/1 = as reference maps. Maps calculated at 6^ = 10~^ should not deviate from the 

reference maps by more than 10^^ in order to declare our error estimate successful. On the other 
hand, we do not want the deviation to be too small either, because this would mean that we are 
making more calculations than required for matching our target accuracy. In this subsection we 
are not considering limb darkening because we want to focus on the accuracy of the single contour 
calculation. 

As a first example, we shall consider a binary lens with mass ratio q = 0.1 and separation 
d = 0.95 (intermediate caustic topology). Fig. [2^ shows the reference magnification map obtained 
with a source radius p^, = 0.01 and a step size Ay = 0.0025 on the source plane in both directions. 
The caustic is very clearly visible, with a spike in the y = position, where we get the maximum 
magnification. 




(55) 



6. Numerical examples 



6.1. Testing the error estimate 
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As a second example, we consider a planetary lens with q = 0.001 and d = 0.95 (resonant 
caustic topology). We choose the resonant caustic topology since it corresponds to the maximum 
caustic extension, allowing finer and more stringent tests. The reference map for p^, = 0.001 and 
Ay = 0.001 is shown in Fig. [2}3. The central spike is much higher in this case, because the source 
radius is much smaller (which is necessary for a better probing of the caustic). 

Now, let us come to the first test. In Fig. [3^ we show the magnification difference map 
between a map calculated at 5fi = W'^ and the reference map shown in Fig. [2^. We can note that 
deviations tend to be spatially correlated far from the caustic, while they are very noisy at caustic 
crossings. However, as it is evident from Fig. [3t, the target accuracy is fully achieved by all points 
in the magnification map. There are just six points with 5fj, > 10"'^, with the maximum error being 
S/j, = 0.012. We can consider this number of points with slightly exceeding error acceptable. We can 
also note that higher magnification points in the map tend to have smaller errors, with deviations 
staying one order of magnitude less than the required accuracy. Without the parabolic correction 
and the exit prescription described at the end of Section HI the discrepancy between the errors of 
low and high magnification points would be much higher, so we consider this as a good result of 
our error estimate strategy. We can barely see something like a damped oscillatory behavior of the 
plot as a function of the magnification. Finally, in Fig. [3)3 we plot the number of sampling points 
versus the magnification. The number of sampling points needed for matching a fixed accuracy 5fi 
grows almost linearly with magnification, which is what we expect in a Green's theorem approach. 

Reducing the source size from = 0.01 to = 0.001 has a slightly beneficial effects on the 
errors, which however stay at the correct order of magnitude, as we can see from Fig. HI 

Coming to the planetary lens, the error map is shown in Fig. [5^, where errors appear much 
more scattered and less concentrated on the caustic. The errors stay at the correct order of 
magnitude and everything seems to be very stable with respect to the mass ratio. 

Finally, we come back to the mass ratio q = 0.1 and p^, = 0.01 and try a map with higher target 
accuracy 5fj, = 10^^. Fig. [6] shows that our sampling strategy and our error estimate performs 
in a very successful way at any target accuracy, with all deviations having the correct order of 
magnitude. 

Since the plots of Fig. [3] and [6] only differ for the target accuracy, it is interesting to see how 
many sampling points we need to add to go from Sfi = 10^^ to 6fi = 10^^. In Fig. [7] we plot the 
ratio between the number of sampling points with 6fi = 10~^ and the number n2 of sampling 
points with 5fj, = 10~^ versus the magnification. Firstly, we see that the ratio does not depend on 
the magnification. Secondly, we can see that the number of sampling points is roughly doubled in 
order to increase the accuracy by a factor of 10. More precisely, the average factor < n^/n2 >= 1.77 
in our maps. This can be understood analytically as follows. Considering that the residual error of 
the parabolic correction goes as but the number of sampling points in the interval [0, 27r] goes 
as A0~^, the accuracy in the magnification goes as n~'*. If n is doubled, the accuracy is improved 
by a factor 16 (it would be just = 4 without the parabolic correction). With < n'^/n2 >= 1-77 
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we have < 6fi2/SfJ'3 >=< ?^3/?^-2 > = 9.81, which is very close to the ratio of the target accuracies 
of the two maps. 

6.2. Linear vs Parabolic approximation 

In the previous subsection we have demonstrated that our error estimators give us a full control 
of the accuracy of our calculations. Now we can go into more detail and try to evaluate the speed-up 
due to the parabolic correction. 

In order to give an estimate as realistic as possible, we should consider the calculation of typical 
microlensing light curves, i.e we should sum up the time spent for calculating the magnification 
along straight lines in the source plane. A good sample of straight lines is provided by the rows 
of our magnification maps. In fact, each row at constant 1/2 can be considered as a straight source 
trajectory parallel to the yi axis with impact parameter uq = y2- Let us denote the number of 
sampling points on each row by npar{uo). This number is proportional to the total time spent for 
calculating the full row at 7/2 = uq. 

Similarly, we shall denote the number of sampling points in the analogous linear calculation 
(without the parabolic correction) by n;j„(no). In order to compare analogous calculations, the 
target accuracy 6fi must be the same in both cases. As error estimator in the linear case, we use 
the parabolic correction itself, which is already available in our code. 

The ratio nun/npar is thus a measure of the speed-up obtained by the introduction of the 
parabolic correction. Fig. [8] shows this quantity as a function of uq for a target accuracy of 
6n = 10~^ (upper points) and 6fi = 10^^ (lower points). The number of sampling points drops 
by a factor 3.3 in average for Sfi = 10~^ and 6.1 for dfj, = 10~^. The speed-up is even higher for 
central events with uq ~ 0, in which high-magnification points have a considerable weight. Such 
numbers are very encouraging, since a factor 6 may bring the computational time e.g. for a huge 
Markov chain from one weak to a single day. 

6.3. Uniform vs Optimal sampling 

The next innovation proposed in this paper is the optimal sampling driven by a reliable estimate 
of the residual error in each arc between two sampling points. Indeed, by increasing the sampling 
only where really needed (e.g. close to caustic crossings of the source boundary), we expect to save 
a considerable amount of computational time. 

In order to evaluate the speed-up, we adopt the same strategy explained in the previous 
subsection: we sum up the number of sampling points for each point in the magnification map 
at fixed y2, so that we can have a measure of the time needed to calculate a microlensing light 
curve for a source trajectory parallel to the yi axis and with uq = 2/2- We denote the number of 



- 20 - 



sampling points with the optimal sampling strategy by Uopt and the number of sampling points with 
a uniform sampling by riuni- The number of sampling points in the uniform case is determined by 
doubling the initial two-points sampling {0, tt} until the new magnification differs from the previous 
one by less than 5^/2, where 6fi is the target accuracy. Note that this prescription is sometimes 
unsafe, since we can have very small features in the images, which could be completely missed. 
However, for the purpose of a gross estimate of the speed-up, we adopt this prescription for the 
uniform sampling, since at most we are just underestimating the correct speed-up in some points. 

In Fig. [9] we plot the ratio nuni/n-opt as a function of uq for trajectories parallel to the yi axis. 
The binary lensing geometry is that described in Fig. [2^. We can see that the speed-up reaches 20 
for central trajectories and then falls down for larger impact parameters. Indeed, when the source 
is poorly magnified, there is no need for optimal sampling. The average speed-up is 3.3 for a target 
accuracy 6fi = 10~^ and 4.8 for (5/x = 10~^, though high magnification events get a much larger 
benefit from optimal sampling. 

6.4. Limb Darkening 

As a final test of our code, we consider a linear limb darkened source with a = 0.51 and radius 
= 0.01. We have generated a reference magnification map with target accuracy 6fj, = 10~^. We 
do not show it because it looks very similar to Fig. [2^, except for the height of the magnification 
peaks. After that, we have generated a test magnification map with target accuracy 6fi = 10~^. 
In both cases we have used the error estimators introduced in section [5] and the optimal source 
sampling strategy described there. 

The difference between the test and the reference map is shown in Fig. [TOb . We can see 
that errors are kept well under control, in particular at caustic crossings, which represent the most 
crucial tests for limb darkening. Of course, approximating limb darkening by the contour method 
requires a large number of points. In Fig. llOb we have a plot of the number of sampling points 
(adding those in all annuli) vs the magnification. We can see that the number of points grows faster 
than linearly with magnification at large ^. The number of annuli (Fig. [TOh ) grows rapidly from 
1 to 20 at low magnifications and then stays more or less constant up to very high magnifications. 
Finally, Fig. [TOl i shows that the target accuracy has been achieved by all points save for three with 
6n = 0.011. 

The bottom line of our limb darkening code is the slow-down plot shown in Fig. [TTJ We 
compare the number of points required for a calculation of limb darkening by summing up the 
number of sampling points in all annuli. We denote this number by ridark- We compare this number 
with the number of sampling points in the uniform source case {rino-dark)- In panel (a) we show 
the ratio ndark/^no-dark vs the magnification. Furthermore, as in previous subsections, we sum up 
the total number of sampling points along each row at fixed y2 in our magnification map. This is 
an indication of the total time spent for calculating a full microlensing trajectory parallel to the yi 
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axis with uq = y2- The ratio of the number of points including the hmb darkening treatment rifiark 
and the number of points calculated for a uniform brightness disk n^o-dark therefore represents the 
slow-down factor of our code for the inclusion of the limb darkening treatment. From Fig. Illb we 
see that the average slow-down factor is 2.64, but central trajectories may be slowed down to a 
factor of 11. Interestingly, we also have two peaks at impact parameters uq = 0.34 and uq = 0.42, 
which correspond to trajectories including cusp crossings. 

7. Conclusions 

In this paper we have presented four innovations for codes attempting microlensing calculations 
based on contour integration. In this class of codes the contours of the images are obtained by 
inversion of the lens map at the source boundary; the area of the images is then obtained by a 
simple contour integration rather than by a surface integration. 

We have introduced a parabolic third order correction in the evaluation of Green's line integral, 
which leaves a residual of the fifth order in the interval size A^. The speed-up with respect to 
the classical trapezium approximation is around 4 for a target accuracy in the magnification of 
5^i = 10-2. 

We have introduced accurate error estimators, whose reliability has been shown by comparing 
several magnification maps obtained at different levels of target accuracy. 

Thanks to these error estimators, we have proposed a new optimal sampling strategy driven by 
the error estimates in each sampling interval. With respect to a uniform sampling we get a speed-up 
ranging from 3 to 20 at 5^ = 10"^ depending on how large is the magnification experienced by the 
source in its microlensing trajectory. 

Finally, we have faced the problem of limb darkening, which is the hardest obstacle for codes 
based on contour integration. Also in this case we have introduced error estimators and an optimal 
sampling strategy with the aim of minimizing calculation keeping full control of the accuracy. As a 
result, we have a very reliable limb darkening approximation to any desired accuracy. At 5^ = 10"^ 
the slow-down with respect to a uniform brightness disk ranges from 2 to 11. 

Summing up, the speed-up gained by the optimal sampling is sufficient to compensate the 
slow-down due to the migration from uniform brightness disks to realistic limb darkened sources. 
In addition, the parabolic correction guarantees a net speed-up with respect to traditional linear 
Green's theorem codes. Finally, the full control of the errors is an invaluable help in order to avoid 
redundant calculations and concentrate the efforts where it is really needed. Its impact with respect 
to traditional codes with fixed sampling is definitely huge and difficult to quantify. 

All these innovations should put codes based on contour integration in the front line for binary 
and planetary microlensing events modelling. 
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Fig. 1. — Linear limb darkening profile with a = 0.51 (dotted line) compared with a block approx- 
imation with four bins (dashed line) and 16 bins (solid line). 
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Fig. 2. — (a) Magnification map for a binary lens with mass ratio q = 0.1 and separation d = 0.95; 
the source radius is = 0.01; the step is Ayi^2 = 0.0025. (b) Magnification map for a binary lens 
with mass ratio q = 0.001 and separation d = 0.95; the source radius is = 0.001; the step is 
Ayi,2 = 0.001. 
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Fig. 3. — (a) Map of the errors for a binary lens with mass ratio q = 0.1 and separation d = 0.95; 
the source radius is = 0.01; the target accuracy is Sfi = 10~^. (b) Number of points used in 
samphng vs magnification, (c) Matched accuracy vs magnification. 
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Fig. 5. — Same as Fig. [3] with a source radius = 0.001 and a planetary mass ratio q = 0.001. 



- 29 - 



u 0.0005^^ 

o.ooooV 

I 

-0.5 




700 
600 
500 

400 
300 
200 t 
]00 






(b) 



• • • • 




20 



40 60 



80 100 



0.0008 - » 
0.0006 - 
0.0004 - 
0.0002 
0.0000 
-0.0002 



Q 



20 



(C) 



•4'** 



40 60 



80 100 



6. — Same as Fig. [3] with a target accuracy Sfi = 10 
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Fig. 7. — Ratio of the number of sampling points for an accuracy 6^ = 10~^ and the number of 
samphng points with 6^ = 10~^. The ratio is plotted versus the magnification. 
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Fig. 8. — Speed-up due to the parabolic correction for microlensing trajectories parallel to the yi 
axis and uq varying from to 0.5 (the reference map is in Fig. [2^). The upper points are for target 
accuracy 6fj, = 10^^ and the lower points are for dfi = 10^^. 
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Fig. 9. — Speed-up due obtained thanks to optimal sampling for microlensing trajectories parallel 
to the yi axis and uq varying from to 0.5 (the reference map is in Fig. [2^). Empty circles are for 
target accuracy 6^ = 10~^; filled circles are for 5fi = 10~^. 
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Fig. 10. — (a) Map of the errors for a binary lens with mass ratio q = 0.1 and separation d = 0.95; 
the source radius is p* = 0.01 with a Hnear limb darkening a = 0.51; the target accuracy is 
(5/Lt = 10~^. (b) Number of points used in sampling vs magnification, (c) Number of annuli vs 
magnification, (d) Matched accuracy vs magnification. 
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Fig. 11. — (a) Ratio between the number of points used with and without the Umb darkening 
treatment, (b) Slow-down due to the hmb darkening treatment for microlensing trajectories parallel 
to the yi axis and uq varying from to 0.5. 



